###############################################################################
#  File-Name:   DDuell_followMajority_runFile.r   		       		
#  Date:        20/05/18                			
#  Author:      Dominik Duell 					    			      						
###############################################################################
pkgs <- c('ggplot2','dplyr','tidyverse','broom','magrittr','coin','clusrank',
          'summarytools','estimatr','car','conflicted','rsample','Hmisc','grid',
          'gridExtra','ggpubr','ggforce','xtable','lme4','stargazer')
for (pkg in pkgs) library(pkg,character.only=T)
conflict_prefer("filter","dplyr")
conflict_prefer("recode","dplyr")
conflict_prefer("lag","dplyr")
options(warn=-1)
###############################################################################
###############################################################################
# Get and prepare data, reorder factor levels as desired, create subsetted data 
# to shorten code further down
###############################################################################
data <- read.csv('DDuell_followMajority_data.csv') %>%
  mutate(voteWinner = recode(voteWinner,`1`='P wins',`2`='R wins'),
    RWins=ifelse(voteWinner=='R wins',1,0),
    treatment=factor(treatment,levels=c('No appeal','Group appeal','Income appeal',
      'Poor MJ - No appeal','Poor MJ - Group appeal','Poor MJ - Income appeal',
      'Rich MI - Group appeal')),
    groupHeterogeneity=factor(groupHeterogeneity,levels=c('Low heterogeneity',
      'Medium heterogeneity','High heterogeneity')),
    groupHeterogeneity.binary=factor(ifelse(groupHeterogeneity=='Low heterogeneity',
      'Low','Medium or high')),
    majorityGroupRich=factor(factor(ifelse(majorityGroup=='MI','MI',
      ifelse(majorityGroup=='MJ'&rich==1,'MJ, rich','MJ, poor')))))

# Main treatments only: No appeal, group appeal, and income appeal treatment
data_richMajority <- data %>% 
  filter(treatment=='No appeal'|treatment=='Group appeal'|
           treatment=='Income appeal') %>% droplevels()

# Main treatments and majority group MJ only, adding indicator variables for 
# coordination mechanism
data_coordinationMechanism <- data_richMajority %>% 
  filter(majorityGroup=='MJ'&groupHeterogeneity!='Low heterogeneity') %>% 
  group_by(groupTreatment) %>% summarise(m.voteR=mean(voteR,na.rm=T)) %>% 
  left_join(data_richMajority %>% filter(majorityGroup=='MJ'&
    groupHeterogeneity!='Low heterogeneity')) %>% 
  mutate(
    m.voteR.noAppeal=ifelse(treatment=='No appeal',
      mean(data$voteR[data$treatment=='No appeal'&data$majorityGroup=='MJ'&
      data$groupHeterogeneity!='Low heterogeneity']),
      ifelse(treatment=='Group appeal',
             mean(data$voteR[data$treatment=='Group appeal'&
                               data$majorityGroup=='MJ'&
                               data$groupHeterogeneity!='Low heterogeneity']),
             mean(data$voteR[data$treatment=='Income appeal'&
                               data$majorityGroup=='MJ'&
                               data$groupHeterogeneity!='Low heterogeneity']))),
    propensity=factor(ifelse(m.voteR>m.voteR.noAppeal,'Propensity to coordinate on R',
      'Propensity to coordinate on P')),
    propensityR=ifelse(propensity=='Propensity to coordinate on R',1,0),
    propensityVsNoAppeal=factor(ifelse(treatment!='No appeal'&
      propensity=='Propensity to coordinate on R',
      'Propensity to coordinate on R',ifelse(treatment!='No appeal'&
      propensity=='Propensity to coordinate on P',
      'Propensity to coordinate on P','No appeal'))))

# Main treatments only, aggregate variables to society-round observations, add 
# strategy profile variables
data_by_society <- data_richMajority %>%
  group_by(groupTreatment,majorityGroup,period,voteR) %>%
  summarise(count = n()) %>% 
  spread(voteR,count) %>%
  mutate(countVoteR=ifelse(is.na(`1`)==F,`1`,0)) %>%
  select(-c(`0`,`1`)) %>%
  spread(majorityGroup,countVoteR) %>%
  mutate(countVoteRMJ = MJ,
    countVoteRMI = MI,
    countVotePMJ = 3-MJ,
    countVotePMI = 2-MI) %>%
  select(-c(MJ,MI))

strategy_profiles <- c('P wins, all vote P','P wins, MJ or MI split',
                       'R wins, MJ or MI split',
                       'R wins, MJ votes R and MI votes P')

# Main treatments only, merge individual-round and society-round-level data s
# set, add equilibrium indicator variables
data_richMajority_equilibria <- left_join(data_richMajority,data_by_society) %>%
  mutate(
    equilibrium = factor(
      ifelse(countVoteRMI==0&countVoteRMJ==0,'(P,P,P;P,P)',
      ifelse(countVoteRMI==0&countVoteRMJ==3,'(R,R,R;P,P)','No equilibrium')),
      levels=c('No equilibrium','(R,R,R;P,P)','(P,P,P;P,P)')),
    eqVoteWinner = factor(
      ifelse(equilibrium=='(P,P,P;P,P)','P wins, all vote P',
      ifelse(equilibrium=='No equilibrium'&
               voteWinner=='P wins','P wins, MJ or MI split',
      ifelse(equilibrium=='No equilibrium'&
               voteWinner=='R wins','R wins, MJ or MI split',
      ifelse(equilibrium=='(R,R,R;P,P)','R wins, MJ votes R and MI votes P',NA)))),
      levels=strategy_profiles),
    whoGetsI=
      ifelse((voteWinner=='R wins'&countVoteRMJ>countVoteRMI)|(voteWinner=='P wins'&
        countVotePMJ>countVotePMI),'MJ wins I',
      ifelse((voteWinner=='R wins'&countVoteRMJ<countVoteRMI)|(voteWinner=='P wins'&
        countVotePMJ<countVotePMI),'MI wins I',
      'MJ and MI split I'))) 

data_exitSurvey <- read.csv('DDuell_followMajority_exitSurvey.csv') %>% 
  mutate(
    sid=treatment*1000+session*100+subject,
    treatment=recode(factor(treatment),`1`='No appeal',`2`='Group appeal',
      `3`='Income appeal',`4`='Poor MJ - No appeal',`5`='Poor MJ - Group appeal',
      `6`='Poor MJ - Income appeal',`7`='Rich MI - Group appeal'))

###############################################################################
###############################################################################
# 4 Experimental design
###############################################################################
# Table 1:  Number of societies, subjects, and subject-round observations for 
# appeal and group heterogeneity treatment conditions. Every session includes 
# 8 rounds of the low group heterogeneity, 28 of the medium group heterogeneity, 
# and 4 rounds of the high group heterogeneity treatment.
###############################################################################
out <- data_richMajority %>%
  group_by(treatment,groupHeterogeneity) %>%
  summarise(societies=length(unique(groupTreatment)),
    subjects=length(unique(sid)),
    rounds=length(sid)) %>%
  spread(groupHeterogeneity,rounds) %>% ungroup() %>%
  add_row(treatment='Total',societies=sum(.$societies),
          subjects=sum(.$subjects),
          `Low heterogeneity`=sum(.$`Low heterogeneity`),
          `Medium heterogeneity`=sum(.$`Medium heterogeneity`),
          `High heterogeneity`=sum(.$`High heterogeneity`))
out 
print(xtable(out,type='latex'),file='table1.tex')

###############################################################################
# 5 Results
###############################################################################
# 5.1 Equilibrium predictions 
###############################################################################
# Statistics in the text (p.18/19): Relative 
# frequency of equilibria by group heterogeneity
# Tests and society-level bootstrap of difference 
# in frequency of strategy profiles and equilibria
#################################################
results.equ.play <- data_richMajority_equilibria %>% 
  group_by(treatment,groupHeterogeneity,eqVoteWinner) %>%
  summarise(n=n()) %>% mutate(N=max(cumsum(n)),prop=n/N) 

results.equ.play %>% group_by(groupHeterogeneity,eqVoteWinner) %>% 
  summarise(N=sum(N),n=sum(n)) %>%
  mutate(N=ifelse(groupHeterogeneity=='High heterogeneity',760,N), 
      # fix at N=760 for High heterogeneity to account in the computation 
      # for no observations of any profile with R winning for high heterogeneity
         prop=n/N)

data_richMajority_equilibria %>% group_by(groupHeterogeneity,whoGetsI) %>% 
  summarise(n=n()) %>% mutate(N=max(cumsum(n)),prop=n/N)  

# Hypothesis tests
# Difference in frequency tests (society-level)
out <- data_richMajority_equilibria %>% 
  mutate(
    REq=ifelse(equilibrium=='(R,R,R;P,P)',1,0),
    PEq=ifelse(equilibrium=='(P,P,P;P,P)',1,0)
    ) %>% 
  group_by(treatment,groupTreatment,groupHeterogeneity) %>% 
  summarise(REq=mean(REq),PEq=mean(PEq),RWins=mean(RWins)) 
out.noHigh <- out %>% filter(groupHeterogeneity!='High heterogeneity')
out.noMedium <- out %>% filter(groupHeterogeneity!='Medium heterogeneity')

# R-equilibrium
out.noHigh %$% wilcox.test(REq~groupHeterogeneity,paired=T,alternative='g')
out.noMedium %$% wilcox.test(REq~groupHeterogeneity,paired=T,alternative='g')
# Dealing with zero differences and ties, two-sided test
out.noHigh %$% wilcoxsign_test(REq~groupHeterogeneity,distribution='exact') 
out.noMedium %$% wilcoxsign_test(REq~groupHeterogeneity,distribution='exact') 

# P-equilibrium
out.noHigh %$% wilcox.test(PEq~groupHeterogeneity,paired=T,alternative='l')
out.noMedium %$% wilcox.test(PEq~groupHeterogeneity,paired=T,alternative='l')
out.noHigh %$% wilcoxsign_test(PEq~groupHeterogeneity,distribution='exact') 
out.noMedium %$% wilcoxsign_test(PEq~groupHeterogeneity,distribution='exact') 

# R wins election
out.noHigh %$% wilcox.test(RWins~groupHeterogeneity,paired=T,alternative='g')
out.noMedium %$% wilcox.test(RWins~groupHeterogeneity,paired=T,alternative='g')
out.noHigh %$% wilcoxsign_test(RWins~groupHeterogeneity,distribution='exact') 
out.noMedium %$% wilcoxsign_test(RWins~groupHeterogeneity,distribution='exact') 

# Bootstrapping inference 
# Take bootstrap sample from set of societies and re-compute frequency of strategy 
# profiles
D <- data_richMajority %>% nest(-groupTreatment)
set.seed(01010)
boots <- bootstraps(D,times=1000) 
boot.results.society.level <- 
  map(boots$splits, ~as_tibble(.) %>% unnest %>% 
    group_by(groupHeterogeneity,groupTreatment,majorityGroup,period,voteR) %>%
    summarise(count = n()) %>% spread(voteR,count) %>% 
      mutate(countVoteR=ifelse(is.na(`1`)==F,`1`,0)) %>%
    select(-c(`0`,`1`)) %>% spread(majorityGroup,countVoteR) %>% 
      mutate(countVoteRMJ=MJ,countVoteRMI=MI,
        countVotePMJ=3-MJ,countVotePMI=2-MI,
        RWins=ifelse(countVoteRMJ+countVoteRMI>=3,1,0)) %>% 
    select(-c(MJ,MI)) %>% 
    mutate(
      REq=ifelse(countVoteRMI==0&countVoteRMJ==3,1,0),
      PEq=ifelse(countVoteRMI==0&countVoteRMJ==0,1,0)
      ) %>% group_by(groupHeterogeneity) %>%
    summarise(REq=mean(REq,na.rm=T),PEq=mean(PEq,na.rm=T),
              RWins=mean(RWins,na.rm=T)) %>% 
      mutate(diff.REq.y=ifelse(groupHeterogeneity=='Medium heterogeneity',
                          REq-lag(REq),REq-lag(REq,2)),
        diff.PEq.y=ifelse(groupHeterogeneity=='Medium heterogeneity',
                          PEq-lag(PEq),PEq-lag(PEq,2)),
        diff.RWins.y=ifelse(groupHeterogeneity=='Medium heterogeneity',
                            RWins-lag(RWins),RWins-lag(RWins,2)),)) %>% 
  bind_rows(.id = 'boots') %>% filter(!is.na(diff.REq.y)) %>% 
  group_by(groupHeterogeneity) %>% 
    summarise(
      diff.REq=mean(diff.REq.y),REq.lower=quantile(diff.REq.y,0.025),
      REq.upper=quantile(diff.REq.y,0.975),
      diff.PEq=mean(diff.PEq.y),PEq.lower=quantile(diff.PEq.y,0.025),
      PEq.upper=quantile(diff.PEq.y,0.975),
      diff.RWins=mean(diff.RWins.y),RWins.lower=quantile(diff.RWins.y,0.025),
      RWins.upper=quantile(diff.RWins.y,0.975)
    )
boot.results.society.level

# Is the mixed strategy equilibrium being played? 
out <- data_richMajority_equilibria %>% 
  filter(groupHeterogeneity=='Low heterogeneity'&
           equilibrium=='No equilibrium') %>%
  group_by(groupTreatment,sid,income) %>% summarise(voteR=mean(voteR)) %>%
  mutate(income=factor(income))

out %>% group_by(income) %>% summarise(voteR=mean(voteR))

m <- lm_robust(voteR~income-1,data=out,cluster=groupTreatment)
summary(m)
linearHypothesis(m,diag(5),c(0,0,.38,.62,.74),vcov.m=cluster.vcov(m,out$groupTreatment))
linearHypothesis(m,'income27 = 0',vcov.m=cluster.vcov(m,out$groupTreatment))
linearHypothesis(m,'income38 = 0',vcov.m=cluster.vcov(m,out$groupTreatment))
linearHypothesis(m,'income44 = .38',vcov.m=cluster.vcov(m,out$groupTreatment))
linearHypothesis(m,'income62 = .62',vcov.m=cluster.vcov(m,out$groupTreatment))
linearHypothesis(m,'income73 = .74',vcov.m=cluster.vcov(m,out$groupTreatment))

#################################################
# Figure 2: Distribution of relative frequency of 
# strategy profiles by group heterogeneity and 
# appeal treatments
#################################################
pdf('barPlotStacked_shareEquPlay_byTreatmentAndIncDistrAndVoteWinner.pdf', 
    height=4)
results.equ.play %>%
  mutate(eqVoteWinner=recode(eqVoteWinner,
    'P wins, all vote P'='Redistributive cand. P wins, all vote P',
    'P wins, MJ or MI split'='Redistributive cand. P wins, MJ or MI split',
    'R wins, MJ or MI split'='Wealth-preserving cand. R wins, MJ or MI split',
    'R wins, MJ votes R and MI votes P'=
      'Wealth-preserving cand. R wins, MJ votes R and MI votes P')) %>%
  ggplot(aes(y=prop,x=treatment,color=eqVoteWinner,fill=eqVoteWinner)) +
  geom_col() +
  facet_grid(~groupHeterogeneity) +
  scale_color_manual(values=c('blue','lightblue','indianred1','red')) +
  scale_fill_manual(values=c('blue','lightblue','indianred1','red')) +
  guides(fill=guide_legend(nrow=2),color=guide_legend(nrow=2)) +
  labs(y='Cumulative relative frequency',x='') + 
  theme_bw() +
  theme(
    legend.position='bottom',
    legend.title=element_blank(),
    legend.spacing.x = unit(.25, 'cm'),
    axis.text.x = element_text(angle=45,hjust=1)
  )
dev.off()

###############################################################################
# 5.2 Appeal treatment effects
###############################################################################
# Statistics in the text (p.20): Relative frequency
# of equilibria by group heterogeneity and appeal
# treatment. Tests of difference in frequency of 
# equilibria and vote share of R as well as
# subject-level bootstrap of difference in vote
# share of R
#################################################
results.equ.play %>% 
  filter(groupHeterogeneity=='Low heterogeneity'&
           eqVoteWinner=='R wins, MJ votes R and MI votes P')

# Hypothesis tests
# Difference in frequency of equilibria, R winning candidate (society-level)
# Low heterogeneity
out <- data_richMajority_equilibria %>% 
  filter(groupHeterogeneity=='Low heterogeneity') %>%
  mutate(REq=ifelse(equilibrium=='(R,R,R;P,P)',1,0),
    PEq=ifelse(equilibrium=='(P,P,P;P,P)',1,0)) %>% 
  group_by(groupTreatment,groupHeterogeneity,treatment) %>%
  summarise(REq=mean(REq),PEq=mean(PEq),RWins=mean(RWins))
# R-equilibrium
# No vs group appeal (p-value):
out %>% filter(treatment!='Income appeal') %$% 
  wilcox.test(REq~treatment,alternative='l')$p.value
# No vs income appeal (p-value):
out %>% filter(treatment!='Group appeal') %$% 
  wilcox.test(REq~treatment,alternative='l')$p.value

# P-equilibrium
# No vs group appeal (p-value):
out %>% filter(treatment!='Income appeal') %$% 
  wilcox.test(PEq~treatment,alternative='g')$p.value
# No vs income appeal (p-value):
out %>% filter(treatment!='Group appeal') %$% 
  wilcox.test(PEq~treatment,alternative='l')$p.value

# R wins 
# No vs group appeal (p-value):
out %>% filter(treatment!='Income appeal') %$% 
  wilcox.test(RWins~treatment,alternative='l')$p.value
# No vs income appeal (p-value):
out %>% filter(treatment!='Group appeal') %$% 
  wilcox.test(RWins~treatment,alternative='g')$p.value

# Medium and high heterogeneity
out <- data_richMajority_equilibria %>% 
  filter(groupHeterogeneity.binary!='Low') %>%
  mutate(RRRPP=ifelse(equilibrium=='(R,R,R;P,P)',1,0),
    PEq=ifelse(equilibrium=='(P,P,P;P,P)',1,0)) %>% 
  group_by(groupTreatment,groupHeterogeneity.binary,treatment) %>%
  summarise(RRRPP=mean(RRRPP),PEq=mean(PEq),RWins=mean(RWins))

# (R,R,R;P,P)-profile
out %>% group_by(treatment) %>% summarise(RRRPP=mean(RRRPP))
# No vs group appeal (p-value):
out %>% filter(treatment!='Income appeal') %$% 
  wilcox.test(RRRPP~treatment,alternative='l')$p.value
# No vs income appeal (p-value):
out %>% filter(treatment!='Group appeal') %$% 
  wilcox.test(RRRPP~treatment)$p.value

# Difference in vote share of R (subject-level)
# Income appeal, poor vs rich for majority group (difference and p-value)
data_richMajority %>% 
  filter(treatment=='Income appeal'&majorityGroupRich!='MI') %>% 
  group_by(groupHeterogeneity,majorityGroupRich,sid) %>% 
  summarise(voteR=mean(voteR)) %>%
  group_by(groupHeterogeneity) %>% 
  summarise(poor=t.test(voteR~majorityGroupRich)$estimate[1],
            rich=t.test(voteR~majorityGroupRich)$estimate[2],
            p.value=clusWilcox.test(
              voteR,group=majorityGroupRich,cluster=sid)$p.value,
            paired.p.value=pvalue(wilcoxsign_test(voteR~majorityGroupRich,
                                                  distribution='exact'))) %>%
  mutate(diff=rich-poor)

# No vs group appeal and no vs income appeal for minority group (difference 
# and p-value)
data_richMajority %>% 
  filter(groupHeterogeneity!='Low heterogeneity'&majorityGroupRich=='MI'&
           treatment!='Income appeal') %>% 
  group_by(treatment,groupHeterogeneity,sid) %>% 
  summarise(voteR=mean(voteR)) %>% 
  group_by(groupHeterogeneity) %>% 
  summarise(noAppeal=t.test(voteR~treatment)$estimate[1],
            groupAppeal=t.test(voteR~treatment)$estimate[2],
            p.value=wilcox.test(voteR~treatment,alternative='g')$p.value) %>%
  mutate(diff=groupAppeal-noAppeal)

data_richMajority %>% 
  filter(groupHeterogeneity!='Low heterogeneity'&majorityGroupRich=='MI'&
           treatment!='Group appeal') %>% 
  group_by(treatment,groupHeterogeneity,sid) %>% 
  summarise(voteR=mean(voteR)) %>% 
  group_by(groupHeterogeneity) %>% 
  summarise(noAppeal=t.test(voteR~treatment)$estimate[1],
            groupAppeal=t.test(voteR~treatment)$estimate[2],
            p.value=wilcox.test(voteR~treatment)$p.value) %>%
  mutate(diff=groupAppeal-noAppeal)

# Bootstrapping inference
# Take bootstrap sample from set of subjects and re-compute vote share of 
# candidate R
D <- data_richMajority %>% nest(-sid)
set.seed(01010)
boots <- bootstraps(D,times=1000) 
results.vote.share.subject.clustered <- 
  map(boots$splits, ~as_tibble(.) %>% unnest %>% 
    group_by(treatment,groupHeterogeneity,majorityGroupRich) %>% 
    summarise(mean=mean(voteR,na.rm=T))) %>% 
  bind_rows(.id = 'boots') %>% 
  group_by(treatment,groupHeterogeneity,majorityGroupRich) %>% 
  summarise(y=mean(mean),lower=quantile(mean,0.025),upper=quantile(mean,0.975))
results.vote.share.subject.clustered

#################################################
# Figure 3: Vote share of candidate R by group 
# heterogeneity and appeal treatments for majority
# MJ and minority MI.  
#################################################
pdf('barPlot_meanAndCIVoteR_byTreatmentAndMajorityGroupAndIncDistr.pdf', 
    height=4)
results.vote.share.subject.clustered %>%
  ggplot(aes(y=y,ymin=lower,ymax=upper,x=treatment,color=majorityGroupRich,
             fill=majorityGroupRich)) +
  geom_bar(stat='identity',position=position_dodge(width=.8)) +
  geom_errorbar(position=position_dodge(width=.8),color='lightgray',width=.2) +
  facet_grid(~groupHeterogeneity) +
  scale_color_manual(values=c('gray','black','black')) +
  scale_fill_manual(values=c('gray','black','white')) +
  labs(y='Vote share of candidate R',x='') +
  theme_bw() +
  theme(legend.position=c(.92,.82),
    legend.title=element_blank(),
    legend.spacing.x = unit(.1, 'cm'),
    axis.text.x = element_text(angle=45,hjust=1))
dev.off()

# Efficiency (Footnote 27, p.20), difference in payoff no vs group appeal 
# (p-values):
data_richMajority_equilibria %>%
  mutate(
    decision.payoff=ifelse(RWins==1,income,income*.5+25),
    identity.payoff=ifelse(majorityGroup=='MJ'&whoGetsI=='MJ wins I',10,
                           ifelse(whoGetsI=='MJ and MI split I',5,0)),
    payoff=decision.payoff+identity.payoff,
    identity.payoff.R = ifelse(RWins==1,identity.payoff,NA),
    identity.payoff.P = ifelse(RWins==0,identity.payoff,NA),
    payoff.R = ifelse(RWins==1,payoff,NA),
    payoff.P = ifelse(RWins==0,payoff,NA)
  ) %>% filter(treatment!='Income appeal') %>% droplevels() %>%
  group_by(groupHeterogeneity) %>%
  summarise_at(vars(payoff,identity.payoff),
               ~t.test.cluster(.,sid,treatment)[20,1])

###############################################################################
# 5.3 Coordination mechanism 
###############################################################################
# Statistics in the text (p.22-25): Propensity of 
# MJ to coordinate on R or P and test for 
# coordination mechanism 
#################################################
# Percentage of societies with a propensity to R/P (p.22):
data_coordinationMechanism %>% group_by(treatment,propensity) %>%
  summarise(n=n()) %>% mutate(N=max(cumsum(n)),prop=n/N) 

# Vote share of R in societies with a propensity to R vs P by treatment (p.22):
data_coordinationMechanism %>% group_by(treatment) %>%
  summarise(
      P=t.test(voteR~propensity)$estimate[1],
      R=t.test(voteR~propensity)$estimate[2],
      p.value=clusWilcox.test(voteR,group=as.numeric(propensity),
                              cluster=groupTreatment)$p.value) %>%
  mutate(diff=P-R)

# Vote share of R, no vs group appeal by income for MJ with propensity to 
# coordinate on R (p.22):
data_coordinationMechanism %>% group_by(income) %>%
  filter(treatment!='Income appeal'&
           propensityVsNoAppeal!='Propensity to coordinate on P') %>%
  summarise(
    noAppeal=t.test(voteR~treatment)$estimate[1],
    groupAppeal=t.test(voteR~treatment)$estimate[2],
    p.value=clusWilcox.test(voteR,group=as.numeric(treatment),
                            cluster=groupTreatment)$p.value) %>%
  mutate(diff=groupAppeal-noAppeal)

# Vote share of R for poorest member if MJ (p.23):
data_richMajority %>% filter(treatment=='Group appeal'&
                               majorityGroup=='MJ'&rich==0) %>% 
  group_by(groupHeterogeneity) %>% summarise(voteR=mean(voteR))
# Low vs medium group heterogeneity (p-value):
data_richMajority %>% filter(treatment=='Group appeal'&
                               majorityGroup=='MJ'&rich==0) %>% 
  group_by(groupHeterogeneity,sid) %>% summarise(voteR=mean(voteR)) %>%
  filter(groupHeterogeneity!='High heterogeneity') %$% 
  wilcoxsign_test(voteR~groupHeterogeneity,distribution='exact') 
# Low vs high group heterogeneity (p-value):
data_richMajority %>% filter(treatment=='Group appeal'&
                               majorityGroup=='MJ'&rich==0) %>% 
  group_by(groupHeterogeneity,sid) %>% summarise(voteR=mean(voteR)) %>%
  filter(groupHeterogeneity!='High heterogeneity') %$% 
  wilcoxsign_test(voteR~groupHeterogeneity,distribution='exact') 

#################################################
# Figure 4: Vote share of candidate $R$ by round 
# for the majority $MJ$
#################################################
D <- data_coordinationMechanism  %>% nest(-sid)
set.seed(01010)
boots <- bootstraps(D,times=1000) 
results <- map(boots$splits, ~as_tibble(.) %>% unnest %>% 
                 group_by(treatment,propensity,rich,period) %>% 
                 summarise(mean=mean(voteR,na.rm=T))) %>% 
  bind_rows(.id = 'boots') %>% 
  group_by(treatment,propensity,rich,period) %>% 
  summarise(y=mean(mean),lower=quantile(mean,0.025),upper=quantile(mean,0.975))

pdf('pointrange_voteRVsRound_byTreatmentAndPropensity.pdf', 
    height=5)
results %>% ungroup() %>% 
  mutate(rich=recode(rich,`0`='Poor in MJ',`1`='Rich in MJ')) %>%
  ggplot(aes(y=y,ymin=lower,ymax=upper,x=period,color=rich,fill=rich,
             shape=rich)) +
  annotate("rect",xmin=0,xmax=41,ymin=-.05,ymax=.05,alpha=.2,
           color='orange',fill='lightgray') +
  geom_line(show.legend=F,position=position_dodge(width=1)) + 
  geom_pointrange(position=position_dodge(width=1),size=.05,fatten=20) +
  facet_grid(propensity~treatment) + 
  scale_shape_manual(values=c(21,24)) + 
  scale_color_manual(values=c('black','black')) + 
  scale_fill_manual(values=c('black','white')) + 
  scale_x_continuous(limits=c(0,41),breaks=c(1,10,20,30,40),
                     labels=c(1,10,20,30,40)) +
  scale_y_continuous(breaks=seq(0,.8,.2),labels=seq(0,.8,.2)) +
  labs(y='Vote share of candidate R',x='Round') +
  theme_bw() +
  theme(legend.position='bottom',
    legend.justification='left',
    legend.box.margin=margin(-10,0,-0,0),
    legend.title=element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()) 
dev.off()

# Number of MJ with positive/negative time trend (p.25, footnote 28):
tidy(data_coordinationMechanism %>% group_by(groupTreatment,period) %>% 
    summarise(voteR=mean(voteR)) %>% group_by(groupTreatment) %>% 
    do(m=lm_robust(voteR~period,data=.)),m) %>% filter(term=='period') %>%
  mutate(time.trend=ifelse(estimate>0&p.value<.05,'Positive, significant',
    ifelse(estimate>0&p.value>=.05,'Positive',
    ifelse(estimate<0&p.value<.05,'Negative, significant',
    ifelse(estimate<0&p.value>=.05,'Negative',NA))))
  ) %>% left_join(data_coordinationMechanism) %>% 
  group_by(treatment,propensity) %>%
  mutate(N=n()) %>% group_by(treatment,propensity,time.trend) %>% 
  mutate(n=n(),percent=n/N*100) %>%
  summarise(N=mean(N),n=mean(n),percent=mean(percent)) %>% 
  filter(treatment=='Group appeal')

# Vote share of R in Rich MJ-Group appeal treatment (p.24):
data %>% filter(treatment=='Rich MI - Group appeal') %>% 
  summarise(voteR=mean(voteR))

# Vote share of R for MJ low vs medium/high heterogeneity (p.24):
data_richMajority %>% 
  filter(treatment=='Group appeal'&majorityGroup=='MJ') %>%
  summarise(
    low=t.test(voteR~groupHeterogeneity.binary)$estimate[1],
    mediumHigh=t.test(voteR~groupHeterogeneity.binary)$estimate[2],
    p.value=clusWilcox.test(voteR,group=as.numeric(groupHeterogeneity.binary),
                            cluster=groupTreatment)$p.value) %>%
  mutate(diff=mediumHigh-low)

# Exitsurvey: How many subjects report appeals did not only matter for them but 
# also others (p.25, footnote 29):
data_coordinationMechanism %>% 
  left_join(data_exitSurvey %>% filter(treatment=='No appeal'|
  treatment=='Group appeal'|treatment=='Income appeal') %>% 
    droplevels(),by=c('sid','treatment')) %>%
  filter(treatment=='Group appeal') %>% 
  mutate(appealMatters=
    ifelse(appealMattersToYou==1&appealMattersToOthers==1,
           'Appeal matters to me and others',
    ifelse(appealMattersToYou==1&appealMattersToOthers==0,
           'Appeal matters to me but not others',
    ifelse(appealMattersToYou==0&appealMattersToOthers==1,
           'Appeal matters to others but not me',
    'Appeals do not matter')))) %>% group_by(propensity,appealMatters) %>% 
  summarise(n=n()) %>% mutate(N=max(cumsum(n)),prop=n/N) %>% 
  filter(appealMatters=='Appeal matters to me and others')

# Vote share of R for rich MJ low vs medium/high heterogeneity (p.25):
data_richMajority %>% 
  filter(treatment=='Group appeal'&majorityGroup=='MJ'&rich==1) %>%
  summarise(
    low=t.test(voteR~groupHeterogeneity.binary)$estimate[1],
    mediumHigh=t.test(voteR~groupHeterogeneity.binary)$estimate[2],
    p.value=clusWilcox.test(voteR,group=as.numeric(groupHeterogeneity.binary),
                            cluster=groupTreatment)$p.value) %>%
  mutate(diff=low-mediumHigh)

###############################################################################
# 5.4 Robustness of experimental results
###############################################################################
out <- data_richMajority_equilibria %>% 
  filter(groupHeterogeneity!='High heterogeneity') %>%
  mutate(thirds=factor(ifelse(period<=12,'Round 1-12',
      ifelse(period>12&period<=24,'Round 13-24',
      ifelse(period>24&period<=36,'Round 25-36',NA)))),
    REq=ifelse(equilibrium=='(R,R,R;P,P)',1,0)) %>%
  filter(groupHeterogeneity=='Low heterogeneity')

# Relative frequency of (R,R,R;P,P) strategy profile by thirds of rounds 1-36 (p.25):
out %>% group_by(thirds,groupHeterogeneity,eqVoteWinner) %>% summarise(n=n()) %>% 
  mutate(N=max(cumsum(n)),prop=n/N) %>% 
  filter(eqVoteWinner=='R wins, MJ votes R and MI votes P')

# Difference in relative frequency of of (R,R,R;P,P) strategy profile first/second 
# vs last third (p.25, p-value):
out %>% group_by(thirds,groupTreatment) %>% summarise(REq=mean(REq)) %>% 
  ungroup() %>% filter(thirds!='Round 13-24') %>% 
  summarise(
    first=t.test(REq~thirds)$estimate[1],
    last=t.test(REq~thirds)$estimate[2],
    p.value=
      pvalue(wilcoxsign_test(REq~as.numeric(thirds),distribution='exact'))) %>%
  mutate(diff=last-first)
  
out %>% group_by(thirds,groupTreatment) %>% summarise(REq=mean(REq)) %>% 
  ungroup() %>% filter(thirds!='Round 1-12') %>% 
  summarise(
    second=t.test(REq~thirds)$estimate[1],
    last=t.test(REq~thirds)$estimate[2],
    p.value=
      pvalue(wilcoxsign_test(REq~as.numeric(thirds),distribution='exact'))) %>%
  mutate(diff=last-second)

# Vote share of R in no appeal and group appeal treatments (main treatments) vs 
# poor MJ-no appeal and poor MJ-group appeal treatments (supplemental 
# treatments, p.26):
data %>% filter((treatment=='No appeal'|treatment=='Group appeal'|
  treatment=='Poor MJ - No appeal'|treatment=='Poor MJ - Group appeal')) %>% 
  group_by(treatment) %>% summarise(voteR=mean(voteR)) 

# Difference in votes share main vs supplemental treatments (no appeal and group 
# appeal, p.26, p-value):
data %>% mutate(
  reverseDistribution=ifelse(treatment=='Poor MJ - No appeal'|
                               treatment=='Poor MJ - Group appeal',1, 
    ifelse(treatment=='No appeal'|treatment=='Group appeal',0,NA))) %>% 
  filter(!is.na(reverseDistribution)) %>% summarise(
    mainTreatments=t.test(voteR~reverseDistribution)$estimate[1],
    supplTreatments=t.test(voteR~reverseDistribution)$estimate[2],
    p.value=
      clusWilcox.test(voteR,group=reverseDistribution,cluster=sid)$p.value) %>%
  mutate(diff=supplTreatments-mainTreatments)

###############################################################################
###############################################################################
# Online appendix
###############################################################################
# B Experimental design appendix
###############################################################################
# B.3 Treatments
###############################################################################
# Table B.2: Summary of all treatment conditions and treatment statistics.
###############################################################################
out <- data %>% group_by(treatment) %>%
  summarise(
    N.societies=length(unique(groupTreatment)),
    N.subjects=length(unique(sid)),
    N=length(sid),
    N.low.het=length(sid[groupHeterogeneity=='Low heterogeneity']),
    N.medium.het=length(sid[groupHeterogeneity=='Medium heterogeneity']),
    N.high.het=length(sid[groupHeterogeneity=='High heterogeneity']))

out
print(xtable(out,type='latex'),file='tableB2.tex')

###############################################################################
# Table B.3: Treatment balance: summary statistics of exit-survey responses
###############################################################################
var.list <- c('age','german','welfare','taxForEducation','taxForWelfare',
              'closeToGroup','groupIDKlee','correctGroupID')
out <- data %>% left_join(data_exitSurvey) %>%
  select(treatment,age,ethnicity,welfare,publicSpending,closeToGroup,klee,
         whichGroupID,groupID,sid) %>% 
  mutate(
    german = ifelse(ethnicity=="Bayern"|ethnicity=="Bayrisch"|
      ethnicity=="Berlin"|ethnicity=="Berlin (egal ob west oder ost)"|
        ethnicity=="Brandenburg"|
      ethnicity=="Brandenburg-Berlin-Region, aber auch Os"|ethnicity=="Deusch"|
        ethnicity=="Deutsch"|
      ethnicity=="Deutsch Griechisch (Europ„isch)"|ethnicity=="Deutsche"|
        ethnicity=="Deutscher"|
      ethnicity=="Deutschland"|ethnicity=="Deutschland, christlich, weiá"|
      ethnicity=="Geburtsort. In meinem Fall: Leipzig"|ethnicity==" Hessen"|
        ethnicity=="Leverkusen"|
      ethnicity=="Mein Elternhaus, deutsch, Tendenz europ Muddastadt Berlin"|
        ethnicity=="Niedersachsen"|
      ethnicity=="Norddeutschland"|ethnicity=="berlin"|ethnicity=="berlin lichtenrade"|
      ethnicity=="berliner"|ethnicity=="de"|ethnicity=="deutsch"|
      ethnicity=="deutsch und sorbisch als Sprache"|ethnicity=="deutsch, europ„isch"|
      ethnicity=="deutschland"|ethnicity=="europ„isch"|ethnicity=="gesamtdeutsch",1,
      ifelse(ethnicity=='',NA,0)),
    taxForEducation=ifelse(publicSpending=='Bildung',1,ifelse(publicSpending!='',0,NA)),
    taxForWelfare=ifelse(publicSpending=='Soziales',1,ifelse(publicSpending!='',0,NA)),
    groupIDKlee=ifelse(groupID==klee,1,ifelse(!is.na(groupID),0,NA)),
    correctGroupID=ifelse(whichGroupID==groupID,1,ifelse(!is.na(whichGroupID),0,NA))
    ) %>% group_by(treatment,sid) %>% summarise_at(vars(var.list),mean) %>% ungroup()
  
# Balance variables (online appendix, p.4, p-value):
out %>% filter(treatment=='No appeal'|treatment=='Group appeal') %>% 
 summarise_at(vars(var.list[-length(var.list)]),
              funs(wilcox.test(.~treatment)$p.value))
out %>% filter(treatment=='No appeal'|treatment=='Income appeal') %>% 
  summarise_at(vars(var.list[-length(var.list)]),
               funs(wilcox.test(.~treatment)$p.value))
out %>% filter(treatment=='No appeal'|treatment=='Poor MJ - No appeal') %>% 
  summarise_at(vars(var.list[-length(var.list)]),
               funs(wilcox.test(.~treatment)$p.value))
out %>% filter(treatment=='No appeal'|treatment=='Poor MJ - Group appeal') %>% 
  summarise_at(vars(var.list[-length(var.list)]),
               funs(wilcox.test(.~treatment)$p.value))
out %>% filter(treatment=='No appeal'|treatment=='Poor MJ - Income appeal') %>% 
  summarise_at(vars(var.list[-length(var.list)]),
               funs(wilcox.test(.~treatment)$p.value))
out %>% filter(treatment=='No appeal'|treatment=='Rich MI - Group appeal') %>% 
  summarise_at(vars(var.list[-length(var.list)]),
               funs(wilcox.test(.~treatment)$p.value))

# Table
count.N <- function(x) sum(!is.na(x))
out.1 <- out %>% group_by(treatment) %>% 
  summarise_at(vars(var.list),list(m=mean,sd=sd,min=min,max=max),na.rm=T) %>% 
  left_join(out %>% group_by(treatment) %>% 
              summarise_at(vars(var.list), funs(N=count.N))) %>% 
  mutate(n=ifelse(treatment!='Rich MI - Group appeal',40,30)) %>%
  pivot_longer(-treatment,names_to=c('variable','.value'),
        names_pattern='(.+)_(.+)') %>%
  filter(!is.na(variable)) %>% select(treatment,variable,N,m,sd,min,max) %>% 
  group_split(treatment)

out.1
print(xtable(out.1[[1]],type='latex',caption='No appeal'),
      file='tableB3_1.tex')
print(xtable(out.1[[2]],type='latex',caption='Group appeal'),
      file='tableB3_2.tex')
print(xtable(out.1[[3]],type='latex',caption='Income appeal'),
      file='tableB3_3.tex')
print(xtable(out.1[[4]],type='latex',caption='Poor MJ – No appeal'),
      file='tableB3_4.tex')
print(xtable(out.1[[5]],type='latex',caption='Poor MJ – Group appeal'),
      file='tableB3_5.tex')
print(xtable(out.1[[6]],type='latex',caption='Poor MJ – Income appeal'),
      file='tableB3_6.tex')
print(xtable(out.1[[7]],type='latex',caption='Rich MI - Group appeal'),
      file='tableB3_7.tex')

###############################################################################
# C Statistical appendix
###############################################################################
# C.1 Summary statistics
###############################################################################
# Table C.6:  Relative frequency of strategy profiles by group heterogeneity 
# and appeal treatments
###############################################################################
out <- results.equ.play %>% spread(groupHeterogeneity,prop) %>% 
  group_by(eqVoteWinner,treatment) %>%
  summarise(Low=mean(`Low heterogeneity`,na.rm=T),
      Medium=mean(`Medium heterogeneity`,na.rm=T),
      High=mean(`High heterogeneity`,na.rm=T)) %>% 
  rename(variable=eqVoteWinner) %>%
  mutate(High=ifelse(is.na(High),0,High)) %>% 
  mutate_at(vars(Low:High),funs(round(.,2))) %>%
  arrange(variable,treatment)

out
print(xtable(out,type='latex'),file='tableC6.tex')

###############################################################################
# Table C.7:  Summary statistics of main variables by income and appeal 
# treatments. Statistics are pooled across all levels of group heterogeneity, 
# subjects, and rounds within one treatment
###############################################################################
out <- rbind(
  data %>% mutate(incomeLevel='All',variable='Vote R') %>% 
    group_by(treatment,variable,incomeLevel) %>% 
    summarise_at(vars(voteR),list(m=mean,sd=sd)) %>% 
    pivot_wider(names_from=treatment,values_from=c(m,sd)),
  data %>% mutate(
    incomeLevel=factor(ifelse(income<30,'Very poor',
      ifelse(income>30&income<50,'Moderately poor',
      ifelse(income>50&income<70,'Moderately rich','Very rich'))),
      levels=c('Very poor','Moderately poor','Moderately rich','Very rich')),
    variable='Vote R') %>%
    group_by(treatment,variable,incomeLevel) %>% 
    summarise_at(vars(voteR),list(m=mean,sd=sd)) %>%
    pivot_wider(names_from=treatment,values_from=c(m,sd)),
  data %>% mutate(incomeLevel='All',variable='R wins election') %>% 
    group_by(treatment,variable,incomeLevel) %>% 
    summarise_at(vars(RWins),list(m=mean,sd=sd)) %>% 
    pivot_wider(names_from=treatment,values_from=c(m,sd)),
  data %>% mutate(incomeLevel='All',variable='income') %>% 
    group_by(treatment,variable,incomeLevel) %>% 
    summarise_at(vars(income),list(m=mean,sd=sd)) %>% 
    pivot_wider(names_from=treatment,values_from=c(m,sd)),
  data %>% mutate(incomeLevel='All',variable='Number of observations') %>% 
    group_by(treatment,variable,incomeLevel) %>%
    summarise(m=length(sid),sd=length(sid)) %>% 
    pivot_wider(names_from=treatment,values_from=c(m,sd)),
  data %>% mutate(incomeLevel='All',variable='Number of subjects') %>% 
    group_by(treatment,variable,incomeLevel) %>%
    summarise(m=length(unique(sid)),sd=length(unique(sid))) %>% 
    pivot_wider(names_from=treatment,values_from=c(m,sd))) %>%
select(variable,incomeLevel,`m_No appeal`,`sd_No appeal`,`m_Group appeal`,
       `sd_Group appeal`,`m_Income appeal`,`sd_Income appeal`,
       `m_Poor MJ - No appeal`,`sd_Poor MJ - No appeal`,
       `m_Poor MJ - Group appeal`,`sd_Poor MJ - Group appeal`,
       `m_Poor MJ - Income appeal`,`sd_Poor MJ - Income appeal`,
       `m_Rich MI - Group appeal`,`sd_Rich MI - Group appeal`)

out
print(xtable(out,type='latex'),file='tableC7.tex')

out %>% filter(variable=='R wins election') %>% 
  select(variable,`m_No appeal`,`sd_No appeal`)
out %>% filter(variable=='Vote R'&incomeLevel=='Moderately poor') %>% 
  select(variable,incomeLevel,`m_Poor MJ - Group appeal`,
         `sd_Poor MJ - Group appeal`)

###############################################################################
# C.2 Additional statistical analysis
###############################################################################
# Table C.8: Multi-level random effects regression of indicator for strategy
# profile (R,R,R;P,P) being played and of indicator for strategy profile 3
# (P,P,P;P,P), P-equilibrium, being played on group heterogeneity treatment, 
# appeal treatment, interaction of those treatments, and round of play 
# including random intercepts for societies.
###############################################################################
# Equilibrium play + relative frequency of vote winner over round
out <- data_richMajority_equilibria %>% mutate(
    RRRPPplayed = ifelse(equilibrium=='(R,R,R;P,P)',1,0),
    PPPPPplayed = ifelse(equilibrium=='(P,P,P;P,P)',1,0))
re1 <- lmer(RRRPPplayed~groupHeterogeneity*treatment+period+
              (1|groupTreatment),data=out)
re2 <- lmer(PPPPPplayed~groupHeterogeneity*treatment+period+
              (1|groupTreatment),data=out)
summary(re1)
summary(re2)
notShow <- capture.output(
  stargazer(re1,re2,column.labels=c('(R,R,R;P,P)','(P,P,P;P,P)'),
            order=c(1,2,3,4,6,7,8,9,5,10),
  covariate.labels = c('Medium heterogeneity','High heterogeneity',
  'Group appeal','Income appeal','Medium heterogeneity $\times$ Group appeal',
  'High heterogeneity $\times$ Group appeal',
  'Medium heterogeneity $\times$ Income appeal',
  'High heterogeneity $\times$ Income appeal','Round','Constant'),
  out='tableC8.tex'))
                                                                                               
###############################################################################
# Figure C.6: Distribution of relative frequency of strategy profiles by group 
# heterogeneity and appeal treatments in first third (round 1-12), second 
# third (round 13-24), and final third (round 25-36) of the experiment. 
# Observations on the high group heterogeneity treatment (round 37-40) are 
# omitted.
###############################################################################
pdf('barPlotStacked_shareEquPlay_byTreatmentAndIncDistrAndVoteWinnerThirds.pdf',
    height=5)
data_richMajority_equilibria %>% 
  filter(groupHeterogeneity!='High heterogeneity') %>%
  mutate(thirds = ifelse(period<=12,'Round 1-12',
    ifelse(period>12&period<=24,'Round 13-24',
    ifelse(period>24&period<=36,'Round 25-36',NA))),
  eqVoteWinner=recode(eqVoteWinner,
    'P wins, all vote P'='Redistributive cand. P wins, all vote P',
    'P wins, MJ or MI split'='Redistributive cand. P wins, MJ or MI split',
    'R wins, MJ or MI split'='Wealth-preserving cand. R wins, MJ or MI split',
    'R wins, MJ votes R and MI votes P'=
      'Wealth-preserving cand. R wins, MJ votes R and MI votes P')) %>% 
  group_by(treatment,thirds,groupHeterogeneity,eqVoteWinner) %>% 
  summarise(n=n()) %>% 
  mutate(N=max(cumsum(n)),prop=n/N) %>% 
  ggplot(aes(y=prop,x=treatment,color=eqVoteWinner,fill=eqVoteWinner)) +
    geom_col() +
    facet_grid(groupHeterogeneity~thirds) +
    scale_color_manual(values=c('blue','lightblue','indianred1','red')) +
    scale_fill_manual(values=c('blue','lightblue','indianred1','red')) +
    guides(fill=guide_legend(nrow=2),color=guide_legend(nrow=2)) +
    labs(y='Cumulative relative frequency',x='') + 
    theme_bw() +
    theme(legend.position='bottom',
      legend.title=element_blank(),
      legend.spacing.x = unit(.25, 'cm'),
      axis.text.x = element_text(angle=45,hjust=1))
 dev.off()
